Setup

knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")


library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)
library(DHARMa)

source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analys
use_mi = FALSE
shutdown = TRUE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  #brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)
df <- openxlsx::read.xlsx('long.xlsx')
df_original <- df

df_list <- prepare_data(df, use_mi = use_mi)

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

df_double <- df_list[[1]]
df_full <- df_list[[2]]

Descriptives

Sample Statistics

# Income (mean of both partner's report)
merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  # Convert to factor
  merged_income <- factor(
    merged_income, levels = c(1,2,3,4,5,6,70,99), 
    labels = c(
      "up to CHF 2'000.-", 
      "CHF 2'001.- to CHF 4'000.-",
      "CHF 4'001.- to CHF 6'000.-",
      "CHF 6'001.- to CHF 8'000.-",
      "CHF 8'001.- to CHF 10'000.-", 
      "above CHF 10'000.-", 
      "I don't know",
      "Undisclosed"
  )
  )
  
  return(merged_income)
}



df_sample_report <- df_full %>%
  group_by(coupleID) %>%
  arrange(userID) %>%
  # Computing couple level variables
  mutate(
    Household_Income = merge_income(first(pre_income_1), last(pre_income_1)),
    
    reldur = pre_rel_duration_m / 12 + pre_rel_duration_y,
    Relationship_duration = mean(reldur, na.rm = TRUE),
    
    habdur = pre_hab_duration_m / 12 + pre_hab_duration_y,
    Cohabiting_duration = mean(habdur, na.rm = TRUE),
    
    Marital_status = factor(
      case_when(
        all(pre_mat_stat == 1) ~ "Married",
        any(pre_mat_stat == 1) ~ "One Partner Married",
        TRUE ~ "Not Married"
      )
    ),
    
    Have_children = factor(
      (first(pre_child_option) + last(pre_child_option)) > 0, 
      levels = c(FALSE, TRUE), 
      labels = c('Have Children', 'No Children')),
    
    Gender = factor(
      gender, 
      levels = c(1,2,3), 
      labels = c('Male','Female', 'Other')),
    

    Couple_type = as.factor(
      case_when(
        first(Gender) == last(Gender) & first(Gender) == 'Male' ~ 'Same-Sex Couple (Male)',
        first(Gender) == last(Gender) & first(Gender) == 'Female' ~ 'Same-Sex Couple (Female)',
        TRUE ~ 'Mixed-sex Couple'
      )
    )
  ) %>%
  ungroup() %>%
  # Individual level variables
  mutate(
    Age = pre_age,
    Handedness = factor(
      pre_handedness, 
      levels = c(0, 1, 2), 
      c('Right','Left', 'Ambidextrous')),
  Highest_Education = factor(
    pre_education, 
    levels = c(1,2,3,4,5,6,7), 
    labels = c(
      "(still) no school diploma",
      "compulsory education (9 years)",
      "vocational training (apprenticeship)",
      "Matura (university entrance qualification)",
      "Bachelor's degree", 
      "Master's degree",
      "Doctorate degree"
      )
    ),
  BMI = pre_weight / ((pre_height / 100)^2) # to meters
  ) %>%
  select(c(Relationship_duration, Cohabiting_duration, Couple_type, Household_Income, 
            Marital_status, Have_children, 
            Gender, Age, Handedness, Highest_Education, BMI))


sample_table <- report_measures(df_sample_report, ICC = F)
sample_table$n_Obs <- as.numeric(sample_table$n_Obs) / 55
rownames(sample_table) <- NULL

n_couple_vars <- 17
sample_table$n_Obs[1:n_couple_vars] <- sample_table$n_Obs[1:n_couple_vars] / 2

packing_sample <- list(
    "Couple level variables (38 couples)" = 
      c(1, n_couple_vars),
    "Individual level variables (76 individuals)" 
    = c(n_couple_vars+1, nrow(sample_table))
    ) 

df_sample_summary <- print_df(
  sample_table,
  rows_to_pack = packing_sample
)

export_xlsx(df_sample_summary,
            file.path('Output', 'SampleDescription.xlsx'),
            merge_option = 'both',
            rows_to_pack = packing_sample,
            colwidths = c(20,35,7,7,7,7,7,10)
            )
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
df_sample_summary
Variable Level n_Obs percentage_Obs Missing Mean SD Range
Couple level variables (38 couples)
Relationship_duration NA 38 NA 0% 9.23 9.03 0.58-36.00
Cohabiting_duration NA 38 NA 0% 7.53 9.14 0.25-33.00
Couple_type Same-Sex Couple (Male) 1 3% NA NA NA NA
Couple_type Same-Sex Couple (Female) 1 3% NA NA NA NA
Couple_type Mixed-sex Couple 36 95% NA NA NA NA
Household_Income I don’t know 0 0% NA NA NA NA
Household_Income up to CHF 2’000.- 2 5% NA NA NA NA
Household_Income CHF 2’001.- to CHF 4’000.- 3 8% NA NA NA NA
Household_Income CHF 6’001.- to CHF 8’000.- 3 8% NA NA NA NA
Household_Income CHF 4’001.- to CHF 6’000.- 5 13% NA NA NA NA
Household_Income Undisclosed 5 13% NA NA NA NA
Household_Income CHF 8’001.- to CHF 10’000.- 8 21% NA NA NA NA
Household_Income above CHF 10’000.- 12 32% NA NA NA NA
Marital_status Married 13 34% NA NA NA NA
Marital_status Not Married 25 66% NA NA NA NA
Have_children No Children 11 29% NA NA NA NA
Have_children Have Children 27 71% NA NA NA NA
Individual level variables (76 individuals)
Gender Other 0 0% NA NA NA NA
Gender Male 38 50% NA NA NA NA
Gender Female 38 50% NA NA NA NA
Age NA 76 NA 0% 34.01 10.96 19.00-60.00
Handedness Ambidextrous 1 1% NA NA NA NA
Handedness Left 11 14% NA NA NA NA
Handedness Right 64 84% NA NA NA NA
Highest_Education (still) no school diploma 0 0% NA NA NA NA
Highest_Education compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education Doctorate degree 0 0% NA NA NA NA
Highest_Education vocational training (apprenticeship) 8 11% NA NA NA NA
Highest_Education Master’s degree 21 28% NA NA NA NA
Highest_Education Bachelor’s degree 23 30% NA NA NA NA
Highest_Education Matura (university entrance qualification) 24 32% NA NA NA NA
BMI NA 76 NA 0% 24.94 4.08 16.37-33.95

Measures

Focal Variables

main_constructs <- c("persuasion", "pressure","pushing", 
                     "pa_sub", "pa_obj", "aff", "reactance"
                     )

main_descriptives <- report_measures(
  data = df_full, 
  measures = main_constructs,
  ICC = TRUE, 
  cluster_var = df_full$userID)

openxlsx::write.xlsx(
  main_descriptives, 
  file.path('Output', 'DescriptivesMain.xlsx')
  )

print_df(main_descriptives)
Variable n_Obs Missing Mean SD Range ICC
persuasion 4180 6% 0.42 1.08 0.00- 5.00 0.23
pressure 4180 6% 0.12 0.62 0.00- 5.00 0.58
pushing 4180 6% 0.16 0.66 0.00- 5.00 0.11
pa_sub 4180 6% 30.40 54.78 0.00-720.00 0.16
pa_obj 4180 11% 144.41 117.81 5.75-971.25 0.18
aff 4180 6% 4.83 1.14 1.00- 6.00 0.41
reactance 4180 81% 0.79 1.32 0.00- 5.00 0.47

All Variables

all_constructs <- c(
  main_constructs,
  "day",
  "weartime",
  "isWeekend",
  "plan",
  "studyGroup",
  "support",
  "got_JITAI",
  "skilled_support"
)


all_descriptives <- report_measures(df_full, all_constructs, ICC = F)

openxlsx::write.xlsx(
  all_descriptives, 
  file.path('Output', 'DescriptivesAll.xlsx')
)

print_df(all_descriptives)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
persuasion NA 4180 NA 6% 0.42 1.08 0.00- 5.00
pressure NA 4180 NA 6% 0.12 0.62 0.00- 5.00
pushing NA 4180 NA 6% 0.16 0.66 0.00- 5.00
pa_sub NA 4180 NA 6% 30.40 54.78 0.00- 720.00
pa_obj NA 4180 NA 11% 144.41 117.81 5.75- 971.25
aff NA 4180 NA 6% 4.83 1.14 1.00- 6.00
reactance NA 4180 NA 81% 0.79 1.32 0.00- 5.00
day NA 4180 NA 0% 0.50 0.29 0.00- 1.00
weartime NA 4180 NA 1% 1057.42 384.29 0.00-1440.00
isWeekend Weekend 1216 29% NA NA NA NA
isWeekend Weekday 2964 71% NA NA NA NA
plan missing 238 6% NA NA NA NA
plan Plan 1860 44% NA NA NA NA
plan No plan 2082 50% NA NA NA NA
studyGroup last 3 weeks interventions 1320 32% NA NA NA NA
studyGroup Allways inerventions 1430 34% NA NA NA NA
studyGroup First 3 weeks interventions 1430 34% NA NA NA NA
support NA 4180 NA 6% 0.91 1.52 0.00- 5.00
got_JITAI JITAI received 585 14% NA NA NA NA
got_JITAI No JITAI 3595 86% NA NA NA NA
skilled_support Days before Intervention 1036 25% NA NA NA NA
skilled_support Days after Intervention 3144 75% NA NA NA NA

Correlations

cors <- wbCorr(df_full[,c(main_constructs)], df_full$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(
  main_cors, 
  file.path('Output', 'Correlations.xlsx')
)

print_df(main_cors, width = '7em')
persuasion pressure pushing pa_sub pa_obj aff reactance
persuasion [0.20] 0.21*** 0.40*** 0.17*** 0.07*** 0.00 -0.05
pressure 0.30 [0.55] 0.28*** -0.03 -0.04* -0.01 0.26***
pushing 0.63*** 0.40* [0.08] 0.14*** 0.05** 0.01 0.07*
pa_sub 0.27 -0.10 0.24 [0.15] 0.31*** 0.20*** -0.04
pa_obj 0.13 -0.08 0.27 0.51*** [0.13] 0.19*** 0.06
aff 0.28 -0.07 0.29 0.52*** 0.22 [0.28] -0.01
reactance 0.18 0.23 0.11 0.06 0.38* -0.12 [0.42]

Within-person correlations are above the diagonal and between-person correlations are below the diagonal.
On the diagonal are intraclass correlations (ICCs)

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,28)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,28+6)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)




prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
  #brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_sub")
)
plot(pa_sub, ask = FALSE)

plot(pp_check(pa_sub, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_sub))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check_transformed(pa_sub, transform = log1p)

loo(pa_sub)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12060.6 177.6
## p_loo        30.2   2.8
## looic     24121.2 355.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.   
## (-Inf, 0.7]   (good)     3732  99.9%  
##    (0.7, 1]   (bad)         3   0.1%  
##    (1, Inf)   (very bad)    1   0.0%  
##             Min. ESS
## (-Inf, 0.7] 662     
##    (0.7, 1] <NA>    
##    (1, Inf) <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 27.35* 21.03 35.95 1.001 3902.26 6229.26
Within-Person Effects
Daily persuasion experienced 1.20* 1.07 1.35 1.000 10192.93 7278.69
Daily persuasion utilized (partner’s view) 1.19* 1.06 1.34 1.000 12159.42 8804.61
Daily pressure experienced 0.94 0.71 1.30 1.000 9685.08 7640.07
Daily pressure utilized (partner’s view) 1.16 0.88 1.59 1.000 11249.31 7539.42
Daily pushing experienced 1.13 0.92 1.41 1.000 8816.39 8383.16
Daily pushing utilized (partner’s view) 1.13 0.95 1.36 1.000 12547.79 9213.16
Day 0.79 0.58 1.10 1.000 12431.32 9455.82
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.61 0.83 3.13 1.002 3172.48 5991.41
Mean persuasion utilized (partner’s view) 1.19 0.62 2.28 1.002 3128.46 5435.36
Mean pressure experienced 0.50 0.22 1.11 1.001 4626.65 6921.35
Mean pressure utilized (partner’s view) 0.43* 0.19 0.96 1.003 4346.50 7133.59
Mean pushing experienced 1.71 0.64 4.56 1.000 4336.92 5727.04
Mean pushing utilized (partner’s view) 2.03 0.76 5.67 1.001 4271.64 6953.82
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.45 0.83 1.00 3236.86 5929.19
sd(Daily persuasion experienced) 0.09 0.00 0.24 1.00 4277.64 4971.03
sd(Daily persuasion utilized (partner’s view)) 0.08 0.00 0.21 1.00 5008.36 5206.23
sd(Daily pressure experienced) 0.17 0.01 0.53 1.00 6480.10 5580.93
sd(Daily pressure utilized (partner’s view)) 0.16 0.01 0.49 1.00 7007.22 4748.59
sd(Daily pushing experienced) 0.28 0.02 0.58 1.00 3266.42 3678.71
sd(Daily pushing utilized (partner’s view)) 0.12 0.00 0.32 1.00 4609.50 3727.29
Additional Parameters
ar[1] 0.03 -0.94 0.94 1.00 10005.25 6639.54
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.00 13938.19 8436.71
sderr 0.05 0.00 0.13 1.00 6131.05 5072.35
sigma NA NA NA NA NA NA

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_obj_log")
)
plot(pa_obj_log, ask = FALSE)

plot(pp_check(pa_obj_log, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_obj_log))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(pa_obj_log, transform = log1p)

loo(pa_obj_log)
## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2813.7  56.5
## p_loo        91.6   4.5
## looic      5627.5 113.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 118.09* 106.05 131.83 1.002 4338.46 6795.25
Within-Person Effects
Daily persuasion experienced 1.03 1.00 1.06 1.001 10995.72 9537.78
Daily persuasion utilized (partner’s view) 1.02 0.99 1.05 1.000 12883.20 8691.54
Daily pressure experienced 0.95 0.89 1.01 1.001 17830.48 9885.92
Daily pressure utilized (partner’s view) 0.98 0.92 1.05 1.000 19278.32 10246.31
Daily pushing experienced 1.03 0.98 1.07 1.000 14025.59 9055.59
Daily pushing utilized (partner’s view) 1.03 0.99 1.07 1.000 21949.21 9623.34
Day 0.96 0.88 1.05 1.000 31435.94 8368.60
Daily weartime 1.00* 1.00 1.00 1.000 11933.20 7885.13
Between-Person Effects
Mean persuasion experienced 1.09 0.81 1.45 1.001 3608.13 6059.46
Mean persuasion utilized (partner’s view) 0.96 0.72 1.29 1.001 3610.97 5890.12
Mean pressure experienced 0.97 0.71 1.34 1.001 5511.94 7801.62
Mean pressure utilized (partner’s view) 0.98 0.72 1.32 1.001 4870.74 7220.18
Mean pushing experienced 1.02 0.66 1.56 1.000 5260.03 7662.48
Mean pushing utilized (partner’s view) 1.29 0.85 1.96 1.000 5271.31 7350.00
Mean weartime 1.00 1.00 1.00 1.000 17022.21 9961.31
Random Effects
sd(Intercept) 0.30 0.23 0.39 1.00 4214.82 7324.07
sd(Daily persuasion experienced) 0.05 0.03 0.08 1.00 7644.01 6901.55
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 0.08 1.00 5351.19 3858.94
sd(Daily pressure experienced) 0.05 0.00 0.14 1.00 6659.89 6610.33
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.11 1.00 7701.49 6409.62
sd(Daily pushing experienced) 0.06 0.00 0.14 1.00 3294.00 5715.05
sd(Daily pushing utilized (partner’s view)) 0.03 0.00 0.08 1.00 6464.76 7907.69
Additional Parameters
ar[1] 0.29 0.26 0.33 1.00 25732.94 7825.38
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.55 0.54 0.57 1.00 22829.74 8539.10

Affect

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "mood")
)
plot(mood, ask = FALSE)

plot(pp_check(mood, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(mood))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(mood, transform = log1p)

loo(pa_sub)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12060.6 177.6
## p_loo        30.2   2.8
## looic     24121.2 355.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.   
## (-Inf, 0.7]   (good)     3732  99.9%  
##    (0.7, 1]   (bad)         3   0.1%  
##    (1, Inf)   (very bad)    1   0.0%  
##             Min. ESS
## (-Inf, 0.7] 662     
##    (0.7, 1] <NA>    
##    (1, Inf) <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.73* 4.52 4.94 1.001 3046.66 4422.02
Within-Person Effects
Daily persuasion experienced 0.00 -0.03 0.04 1.000 17297.90 9842.88
Daily persuasion utilized (partner’s view) 0.02 -0.02 0.05 1.001 14451.84 9886.46
Daily pressure experienced -0.05 -0.16 0.05 1.000 12202.80 8571.79
Daily pressure utilized (partner’s view) -0.04 -0.18 0.09 1.001 11104.70 8944.62
Daily pushing experienced 0.02 -0.04 0.09 1.000 12680.12 9243.11
Daily pushing utilized (partner’s view) 0.06* 0.00 0.12 1.000 13534.43 9546.16
Day 0.22* 0.06 0.38 1.000 24416.16 8358.49
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.32 -0.26 0.88 1.001 2292.03 3556.47
Mean persuasion utilized (partner’s view) 0.23 -0.36 0.80 1.001 2249.94 3776.58
Mean pressure experienced -0.29 -0.89 0.30 1.001 2937.39 4951.55
Mean pressure utilized (partner’s view) -0.31 -0.89 0.28 1.001 2908.62 4194.29
Mean pushing experienced 0.25 -0.55 1.04 1.000 3993.12 6072.68
Mean pushing utilized (partner’s view) 0.39 -0.41 1.18 1.001 3968.95 6710.30
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.59 0.46 0.77 1.00 3439.67 6053.64
sd(Daily persuasion experienced) 0.03 0.00 0.07 1.00 5380.88 6051.41
sd(Daily persuasion utilized (partner’s view)) 0.06 0.01 0.11 1.00 3288.91 3865.59
sd(Daily pressure experienced) 0.11 0.01 0.28 1.00 4275.51 5244.59
sd(Daily pressure utilized (partner’s view)) 0.19 0.03 0.37 1.00 4126.88 4360.12
sd(Daily pushing experienced) 0.09 0.02 0.17 1.00 4973.05 3524.43
sd(Daily pushing utilized (partner’s view)) 0.05 0.00 0.13 1.00 5294.70 5550.35
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 20444.09 8148.47
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 21206.93 8694.13

Reactance

Gaussian

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance")
)
plot(reactance, ask = FALSE)

plot(pp_check(reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance, transform = log1p)

loo(reactance)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1059.5 35.7
## p_loo        75.2  7.5
## looic      2118.9 71.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.8]).
## 
## Pareto k diagnostic values:
##                          Count Pct.   
## (-Inf, 0.7]   (good)     747   98.8%  
##    (0.7, 1]   (bad)        9    1.2%  
##    (1, Inf)   (very bad)   0    0.0%  
##             Min. ESS
## (-Inf, 0.7] 174     
##    (0.7, 1] <NA>    
##    (1, Inf) <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.51* 0.33 0.70 1.000 16417.76 9892.94
Within-Person Effects
Daily persuasion experienced -0.05 -0.11 0.01 1.000 19356.46 9444.65
Daily persuasion utilized (partner’s view) 0.00 -0.06 0.07 1.000 17176.01 9831.82
Daily pressure experienced 0.26* 0.04 0.47 1.000 11824.37 9336.31
Daily pressure utilized (partner’s view) 0.14 -0.06 0.40 1.000 9892.43 7761.10
Daily pushing experienced 0.08* 0.00 0.17 1.000 12338.51 9197.36
Daily pushing utilized (partner’s view) -0.02 -0.10 0.07 1.000 19544.32 8481.36
Day 0.10 -0.17 0.35 1.000 24220.25 9140.35
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.08 -0.29 0.45 1.000 10206.81 9393.41
Mean persuasion utilized (partner’s view) 0.14 -0.25 0.55 1.000 10672.14 9527.96
Mean pressure experienced 0.62* 0.21 1.03 1.000 12564.34 9244.78
Mean pressure utilized (partner’s view) 0.16 -0.28 0.58 1.000 12489.27 9110.28
Mean pushing experienced -0.28 -0.82 0.27 1.000 12347.15 10160.83
Mean pushing utilized (partner’s view) -0.59* -1.16 -0.01 1.000 13098.23 9577.68
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.27 0.15 0.41 1.00 6292.46 7778.91
sd(Daily persuasion experienced) 0.05 0.00 0.13 1.00 3879.53 6215.13
sd(Daily persuasion utilized (partner’s view)) 0.04 0.00 0.13 1.00 7142.65 6976.51
sd(Daily pressure experienced) 0.39 0.22 0.62 1.00 6664.96 9032.62
sd(Daily pressure utilized (partner’s view)) 0.25 0.01 0.63 1.00 3287.32 6408.13
sd(Daily pushing experienced) 0.10 0.01 0.23 1.00 4010.11 5667.82
sd(Daily pushing utilized (partner’s view)) 0.05 0.00 0.15 1.00 8780.84 7897.48
Additional Parameters
ar[1] 0.01 -0.08 0.09 1.00 20581.26 9269.53
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.93 0.88 0.98 1.00 14809.72 8576.95 
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate
## 1 (pressure_self_cw... > 0     0.18
##   Est.Error CI.Lower CI.Upper Evid.Ratio
## 1      0.12    -0.03     0.37      12.61
##   Post.Prob Star
## 1      0.93     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  set_prior("normal(0, 2.5)", class = "b"),
  set_prior("normal(0, 5)", class = "Intercept"),
  set_prior("normal(0, 1)", class = "sd", group = "coupleID", lb = 0),
  set_prior("normal(0, 0.075)", class = "ar", lb = -1, ub = 1),
  set_prior("normal(0.5, 2.0)", class = "sderr", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.95),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance_ordinal")
)
plot(reactance_ordinal, ask = FALSE)

plot(pp_check(reactance_ordinal, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance_ordinal))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance_ordinal, transform = log1p)

loo(reactance_ordinal)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -683.1 32.1
## p_loo       112.0  7.4
## looic      1366.2 64.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.1, 1.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.   
## (-Inf, 0.7]   (good)     753   99.6%  
##    (0.7, 1]   (bad)        3    0.4%  
##    (1, Inf)   (very bad)   0    0.0%  
##             Min. ESS
## (-Inf, 0.7] 143     
##    (0.7, 1] <NA>    
##    (1, Inf) <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA
Intercept[1] 4.29* 2.49 8.21 1.001 2256.01 1951.78
Intercept[2] 9.75* 5.40 21.51 1.002 1435.82 1081.80
Intercept[3] 28.61* 14.76 75.24 1.003 1210.60 778.30
Intercept[4] 133.42* 58.72 466.56 1.004 1094.69 724.17
Intercept[5] 5469.93* 1336.97 45257.45 1.003 1166.93 814.31
Within-Person Effects
Daily persuasion experienced 0.84* 0.70 0.99 1.000 6298.44 5973.36
Daily persuasion utilized (partner’s view) 1.02 0.83 1.24 1.000 10312.32 8638.15
Daily pressure experienced 1.92* 1.20 2.93 1.001 3247.00 2869.57
Daily pressure utilized (partner’s view) 1.25 0.74 2.11 1.000 8335.13 7024.63
Daily pushing experienced 1.17 0.96 1.46 1.000 7133.18 6701.14
Daily pushing utilized (partner’s view) 0.91 0.69 1.19 1.001 9296.13 7955.74
Day 1.48 0.72 3.06 1.000 11603.05 7629.89
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.11 0.38 3.29 1.001 5666.63 6811.76
Mean persuasion utilized (partner’s view) 1.40 0.45 4.79 1.001 5800.97 6413.49
Mean pressure experienced 3.71* 1.20 12.26 1.001 5312.91 5602.60
Mean pressure utilized (partner’s view) 1.20 0.36 3.84 1.001 6106.53 7633.92
Mean pushing experienced 1.18 0.26 5.67 1.000 7661.62 8475.36
Mean pushing utilized (partner’s view) 0.09* 0.01 0.57 1.001 5985.19 7322.95
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.84 0.48 1.31 1.00 2922.60 3336.30
sd(Daily persuasion experienced) 0.18 0.01 0.44 1.00 2882.76 5395.24
sd(Daily persuasion utilized (partner’s view)) 0.22 0.01 0.53 1.00 3661.25 5972.87
sd(Daily pressure experienced) 0.57 0.08 1.15 1.00 2764.88 2568.59
sd(Daily pressure utilized (partner’s view)) 0.45 0.02 1.34 1.00 4034.58 7034.00
sd(Daily pushing experienced) 0.22 0.01 0.52 1.00 3442.46 5129.78
sd(Daily pushing utilized (partner’s view)) 0.19 0.01 0.61 1.00 4713.67 6572.87
Additional Parameters
ar[1] 0.00 -0.15 0.14 1.00 16982.67 8940.39
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.51 0.02 1.65 1.01 613.02 433.90
sigma NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA 
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate
## 1 (pressure_self_cw... > 0     0.49
##   Est.Error CI.Lower CI.Upper Evid.Ratio
## 1      0.26     0.05     0.89      26.84
##   Post.Prob Star
## 1      0.96    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1)
  #brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "is_reactance")
)
plot(is_reactance, ask = FALSE)

plot(pp_check(is_reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(is_reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(is_reactance, transform = log1p)

loo(is_reactance)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -328.4 12.7
## p_loo       260.2 10.7
## looic       656.8 25.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.   
## (-Inf, 0.7]   (good)      35    4.6%  
##    (0.7, 1]   (bad)      572   75.7%  
##    (1, Inf)   (very bad) 149   19.7%  
##             Min. ESS
## (-Inf, 0.7] 396     
##    (0.7, 1] <NA>    
##    (1, Inf) <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.03* 0.00 0.18 1.001 2961.20 5776.43
Within-Person Effects
Daily persuasion experienced 0.55 0.26 1.01 1.000 4130.22 6033.87
Daily persuasion utilized (partner’s view) 1.86 0.76 6.03 1.001 3737.98 5220.40
Daily pressure experienced 8.70* 1.69 60.77 1.000 3710.48 6299.81
Daily pressure utilized (partner’s view) 2.19 0.33 18.42 1.000 6729.25 6651.90
Daily pushing experienced 2.34* 1.07 6.63 1.000 3633.87 4827.24
Daily pushing utilized (partner’s view) 0.72 0.21 2.27 1.001 6444.05 6815.31
Day 3.78 0.38 43.96 1.000 6637.82 7554.77
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 13.42 0.90 291.21 1.000 3985.43 6196.30
Mean persuasion utilized (partner’s view) 5.05 0.31 91.34 1.000 7444.78 7562.53
Mean pressure experienced 132.06* 3.62 5363.94 1.000 9338.38 8352.40
Mean pressure utilized (partner’s view) 6.03 0.15 270.01 1.000 8879.00 9243.83
Mean pushing experienced 4.47 0.10 220.19 1.000 7718.10 8220.27
Mean pushing utilized (partner’s view) 0.06 0.00 3.21 1.000 8490.07 8668.77
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 3.39 1.86 5.30 1.00 2132.54 3149.41
sd(Daily persuasion experienced) 0.58 0.03 1.50 1.00 1978.49 4478.41
sd(Daily persuasion utilized (partner’s view)) 1.50 0.42 2.95 1.00 2329.46 2696.99
sd(Daily pressure experienced) 1.97 0.11 4.47 1.00 1972.68 3600.22
sd(Daily pressure utilized (partner’s view)) 1.35 0.05 3.81 1.00 4817.26 5523.96
sd(Daily pushing experienced) 0.84 0.05 2.02 1.00 2576.31 3993.09
sd(Daily pushing utilized (partner’s view)) 0.75 0.03 2.23 1.00 4299.27 6263.29
Additional Parameters
ar[1] 0.06 -0.16 0.27 1.00 1980.35 3211.48
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 5.90 3.10 9.31 1.00 1442.76 1679.03
sigma NA NA NA NA NA NA 
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate
## 1 (pressure_self_cw... > 0     1.31
##   Est.Error CI.Lower CI.Upper Evid.Ratio
## 1      0.99    -0.23     3.01      11.11
##   Post.Prob Star
## 1      0.92     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random
  model_rows_fixed_final <- model_rows_fixed
  model_rownames_fixed_final <- model_rownames_fixed
  model_rownames_random_final <- model_rownames_random
  rows_to_pack_final <- rows_to_pack
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood,
  reactance,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## [1] "pa_obj_log"
## [1] "mood"
## [1] "reactance"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, 
      "Device-Based MVPA" = 2, 
      "Mood" = 2,
      "Reactance Gaussian" = 2, 
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

summary_all_models
Subjective MVPA
Device-Based MVPA
Mood
Reactance Gaussian
Reactance Dichotome
IRR pa_sub 95% CI pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log b mood 95% CI mood b reactance 95% CI reactance OR is_reactance 95% CI is_reactance
Intercept 27.35* [21.03, 35.95] 118.09* [106.05, 131.83] 4.73* [ 4.52, 4.94] 0.51* [ 0.33, 0.70] 0.03* [0.00, 0.18]
Within-Person Effects
Daily persuasion experienced 1.20* [ 1.07, 1.35] 1.03 [ 1.00, 1.06] 0.00 [-0.03, 0.04] -0.05 [-0.11, 0.01] 0.55 [0.26, 1.01]
Daily persuasion utilized (partner’s view) 1.19* [ 1.06, 1.34] 1.02 [ 0.99, 1.05] 0.02 [-0.02, 0.05] 0.00 [-0.06, 0.07] 1.86 [0.76, 6.03]
Daily pressure experienced 0.94 [ 0.71, 1.30] 0.95 [ 0.89, 1.01] -0.05 [-0.16, 0.05] 0.26* [ 0.04, 0.47] 8.70* [1.69, 60.77]
Daily pressure utilized (partner’s view) 1.16 [ 0.88, 1.59] 0.98 [ 0.92, 1.05] -0.04 [-0.18, 0.09] 0.14 [-0.06, 0.40] 2.19 [0.33, 18.42]
Daily pushing experienced 1.13 [ 0.92, 1.41] 1.03 [ 0.98, 1.07] 0.02 [-0.04, 0.09] 0.08* [ 0.00, 0.17] 2.34* [1.07, 6.63]
Daily pushing utilized (partner’s view) 1.13 [ 0.95, 1.36] 1.03 [ 0.99, 1.07] 0.06* [ 0.00, 0.12] -0.02 [-0.10, 0.07] 0.72 [0.21, 2.27]
Day 0.79 [ 0.58, 1.10] 0.96 [ 0.88, 1.05] 0.22* [ 0.06, 0.38] 0.10 [-0.17, 0.35] 3.78 [0.38, 43.96]
Daily weartime NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.61 [ 0.83, 3.13] 1.09 [ 0.81, 1.45] 0.32 [-0.26, 0.88] 0.08 [-0.29, 0.45] 13.42 [0.90, 291.21]
Mean persuasion utilized (partner’s view) 1.19 [ 0.62, 2.28] 0.96 [ 0.72, 1.29] 0.23 [-0.36, 0.80] 0.14 [-0.25, 0.55] 5.05 [0.31, 91.34]
Mean pressure experienced 0.50 [ 0.22, 1.11] 0.97 [ 0.71, 1.34] -0.29 [-0.89, 0.30] 0.62* [ 0.21, 1.03] 132.06* [3.62, 5363.94]
Mean pressure utilized (partner’s view) 0.43* [ 0.19, 0.96] 0.98 [ 0.72, 1.32] -0.31 [-0.89, 0.28] 0.16 [-0.28, 0.58] 6.03 [0.15, 270.01]
Mean pushing experienced 1.71 [ 0.64, 4.56] 1.02 [ 0.66, 1.56] 0.25 [-0.55, 1.04] -0.28 [-0.82, 0.27] 4.47 [0.10, 220.19]
Mean pushing utilized (partner’s view) 2.03 [ 0.76, 5.67] 1.29 [ 0.85, 1.96] 0.39 [-0.41, 1.18] -0.59* [-1.16, -0.01] 0.06 [0.00, 3.21]
Mean weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 [ 0.45, 0.83] 0.30 [0.23, 0.39] 0.59 [0.46, 0.77] 0.27 [ 0.15, 0.41] 3.39 [ 1.86, 5.30]
sd(Daily persuasion experienced) 0.09 [ 0.00, 0.24] 0.05 [0.03, 0.08] 0.03 [0.00, 0.07] 0.05 [ 0.00, 0.13] 0.58 [ 0.03, 1.50]
sd(Daily persuasion utilized (partner’s view)) 0.08 [ 0.00, 0.21] 0.05 [0.02, 0.08] 0.06 [0.01, 0.11] 0.04 [ 0.00, 0.13] 1.50 [ 0.42, 2.95]
sd(Daily pressure experienced) 0.17 [ 0.01, 0.53] 0.05 [0.00, 0.14] 0.11 [0.01, 0.28] 0.39 [ 0.22, 0.62] 1.97 [ 0.11, 4.47]
sd(Daily pressure utilized (partner’s view)) 0.16 [ 0.01, 0.49] 0.04 [0.00, 0.11] 0.19 [0.03, 0.37] 0.25 [ 0.01, 0.63] 1.35 [ 0.05, 3.81]
sd(Daily pushing experienced) 0.28 [ 0.02, 0.58] 0.06 [0.00, 0.14] 0.09 [0.02, 0.17] 0.10 [ 0.01, 0.23] 0.84 [ 0.05, 2.02]
sd(Daily pushing utilized (partner’s view)) 0.12 [ 0.00, 0.32] 0.03 [0.00, 0.08] 0.05 [0.00, 0.13] 0.05 [ 0.00, 0.15] 0.75 [ 0.03, 2.23]
Additional Parameters
ar[1] 0.03 [-0.94, 0.94] 0.29 [0.26, 0.33] 0.45 [0.42, 0.48] 0.01 [-0.08, 0.09] 0.06 [-0.16, 0.27]
nu NA NA NA NA NA NA NA NA NA NA
shape 0.14 [ 0.13, 0.14] NA NA NA NA NA NA NA NA
sderr 0.05 [ 0.00, 0.13] NA NA NA NA NA NA 5.90 [ 3.10, 9.31]
sigma NA NA 0.55 [0.54, 0.57] 0.87 [0.85, 0.89] 0.93 [ 0.88, 0.98] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::cite_packages()